DNA denaturation bubbles at criticality 
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The equilibrium statistical properties of DNA denaturation bubbles are examined in detail within 
the framework of the Peyrard-Bishop-Dauxois model. Bubble formation in homogeneous DNA is 
found to depend crucially on the presence of nonlinear base-stacking interactions. Small bubbles 
extending over less than 10 base pairs are associated with much larger free energies of formation 
per site than larger bubbles. As the critical temperature is approached, the free energy associated 
with further bubble growth becomes vanishingly small. An analysis of average displacement profiles 
of bubbles of varying sizes at different temperatures reveals almost identical scaled shapes in the 
absence of nonlinear stacking; nonlinear stacking leads to distinct scaled shapes of large and small 
bubbles. 

PACS numbers: 87.10.-e, 87.14.gk, 87.15. Zg 



I. INTRODUCTION 



The nonlinear dynamics and statistical physics of DNA 
denaturation have been widely investigated Recent 
work has attempted to extend mesoscopic scale model- 
ing in order to describe with sufficient accuracy how se- 
quence details determine the statistical and dynamical 
properties of local fluctuations. Such local fluctuations, 
known as "denaturation bubbles" are believed to be in- 
strumental in the initiation of the transcription process 
at physiological temperatures. The possibility of sponta- 
neous, sequence-specific formation of a mid-size bubble 
has been the subject of considerable research interest - 
and some debate - @, H, S S @] • A related - but nonethe- 
less distinct - question which might be relevant to the pro- 
cess of transcription concerns the growth of a bubble to 
much larger sizes. Since this is by definition - at least in 
an asymptotic sense - a scale-free phenomenon, it is best 
addressed at the level of the underlying phase transition; 
moreover, at least its salient features should be evident 
within the context of the homogeneous (polynucleotide) 
chain. 

It should be recalled that nonlinear lattice dynam- 
ics based DNA modeling of the Peyrard-Bishop-Dauxois 
(PBD [7J) type predicts either an (effectively) first-order 
or a (strict) second-order phase transition, depending on 
whether non-linear base-stacking effects are taken into 
account or not Q . In the case of second-order transition, 
it has been determined that domain walls (DW) become 
entropically stable at the critical temperature [{j, [l(| ; re- 
cent numerical evidence from Monte-Carlo simulations 
(Tl| suggests that the average bubble size also becomes 
critical. The first-order transition case is slightly more 
complicated. Entropic effects are not sufficient to enable 
spontaneous DW formation at the critical temperature. 
This appears to rule out DWs as agents of thermal denat- 



uration. Bubbles are natural - and in fact have always 
been - prime suspects for this role. Very recent work 
[l2| has demonstrated that large bubbles may form in 
this case as well, and that their probability distribution 
cannot be described by a simple exponential. Detailed 
data in the vicinity of the critical temperature are not 
available 13[ . It is however known from previous work 
d, 03] that most of the physics of nonlinear base-stacking 
is generated by an effective thermal barrier which modi- 
fies the on-site Morse-like potential. The point at which 
the effects of the thermal barrier become important de- 
fines a natural crossover between different types of be- 
havior. As this note will show in some detail, such a 
crossover is also present in bubble statistics. Bubbles 
which extend over a few sites are entirely non-critical; 
the onset of criticality is reflected only in the statistics of 
large bubbles; the asymptotic properties of the latter are 
such that the free-energy barrier toward bubble growth 
is lowered - approaching zero at the transition temper- 
ature, whereas the average bubble size always remains 
finite. It will be further shown that this dichotomy be- 
tween large and small bubbles is not restricted to the 
statistics; reduced shapes of large and small-size bub- 
bles - which scale uniformly in the absence of nonlinear 
stacking interactions - reveal distinct differences when 
nonlinear stacking is included. 



The paper is organized as follows: Section II includes 
model definitions, notation and general properties of bub- 
bles. Section III presents numerical results on bubble 
statistics based on direct matrix multiplication. Sec- 
tion IV formulates an alternative procedure based on an 
associated eigenvalue problem which provides emphasis 
on asymptotic properties. Section V discusses bubbles 
shapes. The final section includes a brief summary and 
discussion of some key points. 
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II. DEFINITIONS 



Model 



The PBD model assumes a potential energy of the form 



H P = J2{W(y J -uy J ) + V(y J )} 



(1) 
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where yj is a transverse coordinate representing the sep- 
aration of the two bases at the jth site, 
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is an anharmonic elastic term which models the nonlin- 
ear base-stacking interaction and V(y) = (1 — e~ v ) 2 is an 
on-site Morse potential describing the combined effects 
of hydrogen-bonding, stacking and solvent. I will use the 
dimensionless parameter values R = 10.1, b = 0.08 and, 
unless otherwise stated, p = 1. Furthermore, I will as- 
sume that the system is subjected to periodic boundary 
conditions. Thermodynamics is governed by the proper- 
ties of the transfer integral (TI) equation 



dy'K(y,y')My')=^Mv) 
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with K{y,y') = r^WW/T and T the 
dimensionless temperature. In particular, details of a 
possible phase transition depend on the type of singular- 
ity (if any) which the spectral gap Ae = — Tln(Ai/Ao) 
might exhibit near a critical temperature T c . Since the 
spectral gap is equal to the singular part of the ther- 
modynamic free energy per site [9(, a linearly vanishing 
gap as T — > T~ corresponds to a first order transition, a 
quadratically vanishing gap to a second order transition 
etc. 
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where Zn is the full configurational partition function, 
dominated by the highest eigenvalue Ao, 4>o denotes the 
TI eigenstate corresponding to Ao, and K = K/Aq. 

C. Sum rules 

By definition, the sum of all P„'s expresses the proba- 
bility that the site which precedes the bubble has a bound 
base pair, i.e. 



n=0 

Moreover, the sum 
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expresses the fraction of sites with unbound base pairs. 
As a consequence, the average bubble size (including the 
correct weighting factor for bubbles of zero length) is 



1-p 



(8) 



Some general conclusions can already be drawn at this 
level. For a second-order transition, where p oc T c — T 
near T c [8|, oc (T c — 7 1 ) -1 . For a first-order transition, 
where p approaches a constant as T — > T~ , t;b remains 
finite. Note that the average bubble length is in both 
cases much smaller that the correlation length £ = T j Ae 
which respectively diverges quadratically or linearly. 



B. Bubbles 

The nth base pair is assumed to be unbound if y n > y c ; 
it is in a bound state if y < y c ; I choose y c = In 2, 
the inflexion point of the Morse potential. The choice 
is of course somewhat arbitrary, but it should not influ- 
ence fundamental asymptotic results. A bubble of length 
n is a sequence of n successive unbound sites preceded 
and followed, respectively, by a single bound site. It is 
present in the infinite system (assumed subjected to pe- 
riodic boundary conditions) with a probability 
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III. RESULTS 

For bubbles of up to moderately large size, it is possible 
to obtain results by direct matrix multiplication of (O. 
I use a grid of 2989 points and perform the successive 
integrations using a 10th order Bode[l5| routine in the 
interval (-5,205). 



A. Linear base-stacking (p = 0) 

Fig. Q] shows the probability of bubble occurrence (left 
panel) for a variety of temperatures. It is clear that that 
the distribution is far from exponential . For compar- 
ison I show in the inset the probability of occurrence of 
a bound cluster of n sites, which is described by a pure 
exponential. Note that it is possible to fit the P n data, 
at least at the lowest temperatures, by stretched expo- 
nentials, i.e. P n = aexp[—(n/a) b ]. However, the fits are 
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FIG. 1: (Color online) Bubble statistics in the case of vanish- 
ing nonlinear base stacking, p — 0. Left panel: Probability of 
a bubble extending to n sites, for a variety of temperatures. 
Note that the decay becomes weaker as the critical tempera- 
ture is approached. Also shown are, for the two lowest tem- 
peratures, fits to a stretched exponential function. Inset: For 
comparison, I show the probability of a cluster of n succsessive 
sites with bound base pairs (pure exponential). Right panel: 
Ratios of successive probabilities vs. inverse bubble size for a 
range of temperatures. 



neither perfect (they show systematic deviations at both 
ends) nor very instructive (the a obtained decreases as 
the temperature increases, i.e. a — 2.04 at T = 0.85, 
a = 1.07 at T = 0.95 the necessary compensation is 
achieved by a substantial decrease in the stretching ex- 
ponent from b = 0.70 to b — 0.52). A more promising 
approach is to separate out any implicit exponential de- 
pendence by looking at the ratios of successive probabil- 
ities. Thus if one attempts to describe deviations from 
exponential dependence by a power-law, 
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a plot of the ratios 
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vs 1/n should approach a definite limit as n — > oo, from 
which it is in principle possible to read off both the ac- 
tivation free energy Afb = T/a associated with bubble 
growth and - by estimating the asymptotic slope - the 
exponent c [l7} ■ The successive ratios shown in the right 
panel of Fig. [1] confirm this picture; results of the ex- 
trapolation are summarized in Fig. [21 They show that 
it is indeed consistent to represent bubble statistics by 
©, that A/f, ex (T c — T) 2 and that the exponent c varies 
significantly with temperature, approaching a value close 
to 3/2 near T c , 

An appropriate measure of the quality of the numeri- 
cal data is given by the estimated critical temperatures 
(cf. intersections of the dotted lines in Fig. [1] with the 
horizontal axis). Estimates obtained, respectively, from 
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FIG. 2: (Color online) Summary of critical results in the case 
p = 0: as extracted from the numerical TI solution and the 
asymptotic behavior of the curves shown in the right panel 
of Fig. [T] (i) spectral gap Ae (open squares) vs. temper- 
ature, (ii) activation free energy Afb associated with bubble 
growth (open circles); (iii) the square root Ae 1//2 (full squares) 
is known to depend linearly on the temperature; rounding is 
due to the finiteness of the matrix used; (iv) Afl (full cir- 
cles) is found to depend linearly on the temperature; remark- 
ably, it exhibits no rounding; (v) the fraction of bound sites p 
(diamonds), as obtained from (JB}. The dotted lines represent 
linear fits to the data (cf. text for discussion). Inset: the 
exponent c vs temperature. 



the fraction of bound sites p and from A/&, are 1.227(6) 
and 1.229(1). They should be compared with the value 
1.2276(4) obtained via systematic finite-size scaling anal- 
ysis [19j . 



B. Nonlinear base-stacking 



I now proceed to the physically more relevant case of 
nonlinear base-stacking, p = 1. Although this value un- 
derestimates the importance of nonlinear stacking inter- 
actions, it facilitates the present discussion because the 
relevant crossover effects occur in numerically observable 
regions. The important qualitative features remain un- 
changed. In particular, the transition is for all practical 
purposes a first-order one. The TI spectral gap vanishes 
linearly near T c = 0.801 and the fraction p of bound sites 
has an apparent discontinuity at that temperature (cf. 
Fig [3]). In order to avoid data cluttering, the left panel of 
Fig. 0] shows the function P n for two temperatures only, 
T = 0.76 and T — 0.79; note that the higher temperature 
is quite close to the critical temperature. 
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It is possible to obtain rough fits to the full sets of data 
with the functional form © (dotted lines); the parame- 
ter values obtained are a — 8.4, 12.4 and c = 1.42, 1.49 
- where the second value refers to the higher tempera- 
ture. Such fits over the entire data range are of course 
of questionable value if one tries to extract asymptotic 
information. I have included them because the extracted 
parameter values provide a hint of the underlying prob- 
lem: the er's, although much larger than unity are sig- 
nificantly smaller than the correlation length extracted 
from the spectral gap (£ = 19.6, 56.5 for the tempera- 
tures under consideration); more importantly, values of 
c < 2 in conjunction with large a values (indicating the 
onset of criticality) imply - from the known properties of 
the polylogarithm function - a divergent average bubble 
length, which directly contradicts the exact sum rule (cf. 
above). 
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FIG. 3: (Color online) Summary of critical results in the 
case p — 1: (i) spectral gap Ae (diamonds) and (ii) fraction of 
bound sites p (stars) vs. temperature, as obtained from the TI 
numerical solution; (iii) activation free energy Afb associated 
with bubble growth (filled squares) and (iv) exponent c of the 
asymptotic form (I10|) (triangles, inset), as extracted from the 
n — * oo asymptotics of Fig. [4] The dotted lines represent 
linear fits to the data (i) and (ii), both yielding, respectively, 
T c — 0.80 within numerical accuracy. 

The reason behind these difficulties becomes obvious 
if one looks at the ratios P n /P n _i, shown in the right 
panel of FigJU It then becomes clear that small bubbles 
have an entirely different behavior than larger bubbles. If 
one restricts attention to smaller bubbles, e.g. n < 8 for 
T = 0.76 (a temperature quite close to T c ), the apparent 
asymptotic value of the ratio (cf. dotted line) appears 
much smaller; the physics behind this is that it takes a 
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FIG. 4: (Color online) Bubble statistics in the case of nonvan- 
ishing nonlinear base stacking, p = 1. Left panel: Probability 
of a bubble extending to n sites at T = 0.76 and T = 0.79; 
also shown (continuous lines) are fits to the function ©, cf. 
text for a discussion of the fitting parameters. Inset: For com- 
parison I show the probability of a cluster of n bound sites 
(pure exponential). Right panel: Ratios of successive proba- 
bilities vs. inverse bubble size for a range of temperatures. 
Inset: a zoom of the asymptotic region for the two highest 
temperatures; dotted lines represent quadratic extrapolations 
from the last 3 points. 



free energy which is typically higher than T in order to 
generate such a bubble. Note that the low apparent val- 
ues of the exponent c derived from small bubbles (e.g. 
c = 0.75 at T = 0.76) are irrelevant in this context, be- 
cause they are not accompanied by a sequence of appar- 
ent A/b's approaching zero at the critical temperature - 
hence no divergences in the either one of the series © or 
are generated. 

A proper analysis of the asymptotics of large bubbles 
(cf. inset of FigQJ shows that the ratios P n /P n _i lead 
to A/b's with the correct limiting behavior, i.e. linearly 
vanishing at T c (filled squares in Fig. [3]). Moreover, the 
values of c associated with the true asymptotics are also 
much larger (cf. inset in Fig. [3]). At a temperature 
very near T c , an estimate c = 2.9 is obtained, clearly 
consistent with the finite average bubble size demanded 
by a first-order transition (cf. above). 



IV. AN ALTERNATIVE PROCEDURE 

It is possible to use simple linear algebra in order to 
extract some formal properties of the limit of very large 
bubbles n — ► oo. Using the overlap matrix elements 
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it is possible to rewrite (O in the form 
P n =<A\C---C\A> 
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where the matrix product contains n factors, C vv i = 
(A V A'„/A2) 1/2 (<W -B vv ,) and A v = {K v / K^B^. It 
then follows that 



P„ 



(12) 



where {/i Q } are the eigenvalues of the real symmetric 
matrix C and a a = (SA) a , where 5 is the orthogonal 
matrix which diagonalizes C. 

If the spectrum of C is continuous - as the numerical 
computations suggest - then it is possible to rewrite the 
bubble probability in a continuum form 



/Zg / dxTl{x)(i{x)e~ xn ' T 



where the eigenvalues are now labeled as \i a — > /iqc 



(13) 



-x/T 

with a density TZ(x), and |a a | 2 — > /3(x). In the limit 
of very large n only a very narrow range a; < T/n con- 
tributes to the integral. Therefore, if IZ(x) oc a; - ^ and 
f3(x) oc a;'' in the neighborhood of zero, it follows that 
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and, by comparison with ([9]), we conclude that in the 
PBD model the exponent c = r\ — £ + 1 reflects the 
behavior of the density of eigenvalues of C and the /3 
function near the high end of the spectrum. Further- 
more, the highest eigenvalue /iq can be used to extract 
A/ 6 /r = l/<7 = -ln M o. 

Fig. [5] illustrates the behavior of the functions IZ(x) 
and f3(x) in the cases of both linear and nonlinear stack- 
ing. In the first case (left panel), the asymptotics lead 
to a value c = 1.33 for T = 1.20, which is close, but not 
identical with the 1.47 obtained by the extrapolation pro- 
cedure in the previous section. In the nonlinear stacking 
case (right panel), it is seen that the crossover from small 
to larger bubbles has its origins in a strong anomaly of 
the (3{x) function near x = 0.1. Using an effective ex- 
ponent r) = 2.5, and £ = 0.4 would imply c = 3.1, in 
reasonable agreement with the extrapolation estimates 
of the previous section. 



V. BUBBLE SHAPES 

The alternative method described in the previous sec- 
tion, owing to its superior computational efficiency, is 
uniquely suited to deal with the repeated calculations in- 
volved in obtaining full (average) displacement profiles of 
bubbles. A straightforward generalization of the above 
scheme allows the calculation of the conditional average 
displacement of the sth site in a bubble of size n, 



y(s\n) = i-J2' 
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FIG. 5: (Color online) The functions TZ(x) (open symbols, 
left y-axis) and (3{x) (filled symbols, right y-axis). Left panel, 
p = 0: T = 0.85 (squares), T = 1.20 (circles). The density of 
eigenvalues can be fitted with a power law £ = 0.51, regardless 
of temperature; the function f3(x) also exhibits a power-law 
behavior as x — > 0; at the highest temperature a value r\ — 
0.81 is obtained. Right panel: p = 1: T = 0.63 (squares), 
T = 0.79 (circles). The density of eigenvalues exhibits a slight 
anomaly around x = .1; in the region x — * it follows a 
power law £ = 0.41 (slope of dotted line, guide to the eye). 
The function /3(x) exhibits a strong anomaly around x — 0.1 
and does not seem to settle to a pure power law behavior; 
however, the effective slope is quite high (dashed line, guide 
to the eye has a slope of 2.5). 



where D = SDS' 1 and 



dy (t>v{y) y <t>v (y) 



(15) 



Note that the double sum in (fl"4"|) is - apart from a factor 
y s in the integrand - essentially the statistical weight |(SJ| . 

Fig. [6] summarizes the results obtained via this ap- 
proach for bubble shapes in the case of vanishing non- 
linear stacking. The left panel shows that if the reduced 
average displacements (i.e. divided by the maximal dis- 
placement found for each bubble) are plotted against the 
relative site coordinate s/n, the shapes obtained are in- 
dependent of bubble size and/or temperature. Results 
are well fitted by the shape 



y(s\n) = y° n (T)[Ax(l ~ x)Y 



(16) 



where y° (T) can be understood as an average amplitude 
of an n— site bubble and r = 1/2. The dependence of 
bubble amplitudes on size and temperature is shown in 
the right panel. The amplitudes vary with size according 
to a power law with an exponent close to 1/2. The overall 
temperature dependence appears to be roughly linear. 

A note of caution is due at this point. The results de- 
scribed in this paper do not demonstrate the existence 
of a bubble as a well-defined long-lived entity of a given 
fixed size. The interpretation of the shape reported here 
is somewhat more indirect, since it concerns statistical 
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FIG. 6: (Color online) Bubble shapes and sizes in the case 
p — 0. Left panel: Relative average displacements of sites in 
a bubble vs. relative site coordinate. Results from a variety 
of sizes and temperatures collapse on a single curve. Right 
panel: The maximal displacement of an n— site bubble for 
various temperatures. 
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FIG. 7: (Color online) Bubble shapes and sizes in the case 
p = 1. Left panel: Relative average displacements of sites in 
a bubble vs. relative site coordinate. Results from different 
temperatures and sizes collapse on different curves for large 
(open symbols) and small (filled symbols) bubbles. Right 
panel: The maximal displacement of an n— site bubble for var- 
ious temperatures. Note the crossover which occurs around 
n = 10 — 20 at temperatures both near and far from T c . 



average profiles of fluctuating objects. On the other hand 
if a time-dependent solution of a finite extent exists then 
its time-averaged spatial profile, after allowing for correc- 
tions due to interactions with phonons, should look like 
the left panel of Fig. [51 

Long-lived entities with an internal oscillation (discrete 
breathers (DB)), have been reported [20j in the class of 
models under consideration here. Unfortunately, most of 
the work done on DBs concerns objects which are very 
localized in space. It would be interesting to examine 
whether approximate breather-like excitations, perhaps 
with shorter lifetimes and extending over many lattice 
sites, could produce average displacement profiles com- 
patible with (Unj). 

In the case of nonlinear stacking, the results of the 
previous section indicate that small and large bubbles 
behave differently. The derived shapes (Fig. right 
panel) confirm this dichotomy. Small-size bubbles tend 
to have a slightly different shape from large-size bubbles. 
The difference can be expressed in terms of the exponent 
r of the size function (fTOjl . with r = 0.55 for small-size 
and t = 0.82 for large-size bubbles. A similar dichotomy 
occurs when we look at the dependence of bubble am- 
plitude vs. size. There is a crossover behavior from a 
low-exponent to a high-exponent region (cf. the dotted 
and dashed lines with slopes 1/2 and 1 respectively serv- 
ing as guides to the eye). This type of behavior comple- 
ments the findings on the incremental free energy needed 
for bubble growth. Amplitude grows weakly with size 
for small bubbles. After a bubble size threshold (10 — 20 
sites) has been crossed growth becomes a lot easier. 



VI. SUMMARY AND DISCUSSION 

The detailed properties of locally denatured regions 
(bubbles) of homogeneous DNA have been discussed in 
the framework of the PBD model. It has been shown 
that bubble statistics is very sensitive to the presence of 
nonlinear stacking interactions. In summary: 

In the limit of linear base stacking (p = 0) the distri- 
bution of bubble sizes can be described by the product of 
an exponential and a power law (cf. Eq. [9]). The latter 
is characterized by an exponent c which is weakly de- 
pendent on temperature. The exponent's value near the 
critical temperature, close to 3/2, suggests a formal anal- 
ogy with the Poland-Scheraga description of 3-d polymer 
loops in the random walk limit. However, the analogy 
does not rest on microscopic footing; the logarithmic cor- 
rection to the bubble entropy is not related to looping 
and, more importantly, it is not even approximately con- 
stant in temperature; on the other hand, the value of c 
must approach 3/2 at the critical temperature since the 
transition is known to be exactly second order in the case 
P = 0. 

In the presence of nonlinear base stacking (p = 1) the 
analysis of successive probability ratios P n /P n _i reveals 
a more complex behavior of the exponent c. Depending 
on the range of bubble sizes analyzed, different apparent 
values of c are estimated. Small bubbles lead to small 
values of c, larger bubbles suggest larger c values; the lat- 
ter are of course consistent with the apparent first-order 
transition. The threshold seems also to control average 
bubble shapes. Note that it is possible to relate the value 
of the threshold to the parameters of the thermal barrier 
U(y) = (T/2) ln(l + pe^) @, G3 induced by the non- 
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FIG. 8: (Color online) Comparison of Aft, = — Tln/xo ob- 
tained in the context of the associated eigenvalue problem 
with the TI spectral gap Ae in the cases p = 1 (open sym- 
bols) p = (filled symbols); in the latter case the square root 
is plotted. 
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